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An Explicit Scheme for Incorporating Ambipolar Diffusion in a 

Magnetohydrodynamics Code 

Eunwoo Choi 1 , Jongsoo Kim 1,2 , and Paul J. Wiita 3 ' 4 
ABSTRACT 

We describe a method for incorporating ambipolar diffusion in the strong cou- 
pling approximation into a multidimensional magnetohydrodynamics code based 
on the total variation diminishing scheme. Contributions from ambipolar diffu- 
sion terms are included by explicit finite difference operators in a fully unsplit way, 
maintaining second order accuracy. The divergence-free condition of magnetic 
fields is exactly ensured at all times by a flux-interpolated constrained transport 
scheme. The super time stepping method is used to accelerate the timestep in 
high resolution calculations and/or in strong ambipolar diffusion. We perform 
two test problems, the steady-state oblique C-type shocks and the decay of Alfven 
waves, confirming the accuracy and robustness of our numerical approach. Re- 
sults from the simulations of the compressible MHD turbulence with ambipolar 
diffusion show the flexibility of our method as well as its ability to follow complex 
MHD flows in the presence of ambipolar diffusion. These simulations show that 
the dissipation rate of MHD turbulence is strongly affected by the strength of 
ambipolar diffusion. 
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Introduction 



In galactic molecular clouds, ambipolar diffusion, which arises i n partially ionized plas- 



Mouschovias 


1976; 


Shu et al. 


1987) 



(e.g., 



Mestel fc Spitzerl Il956 



lar gas is dense enough that recombination is nearly total, so t hat very low fractions (< 10 



-7^ 



of the gas remains ionized while the rest of the gas is neutral (jMyerslll985l ). This small resid- 
ual ionization is usually attributed to cosmic rays, which can penetrate nearly all clouds. 
Ambipolar diffusion causes the relative drift of ions coupled to the magnetic field and neutrals 
in the molecular cloud cores and so it enable the cloud cores to collapse gravitationally. 

Star formation assisted by ambipolar diffusion has been st udied extensively in the con- 



text of magnetically subcritical or supercritical models (see e.g., 


Mouschovias & Ciolek 


1999 


Desch & Mouschovias 


2001; 


Basu & Ciolek 


2004; 


Tassis & Mouschovias! 


2007 


). In a current 



paradigm of star formation, magnetically supported molecular cloud cores must lose mag- 
netic support through the action of ambipolar diffusion so that star formation can take place. 
Rece nt works have focused on the role of turbulence in the formation of protostellar cores 



^e.g., iNakamura fc Lil 120081 ) . Including the effect of turb ulence on the mechanism of am- 
bipolar diffusion can enha nce the ambipolar diffusion rate (IFatuzzo fc Adamsll2002l ; IZweibel 
20021 ; iHeitsch et al.l 120041 ) , so that the ambipolar diffusion timescale is significantly shorter 
than that estima ted for a similar, but qu iescen t, medium. Usi ng three-dimensional numer- 
ical simulations, lOishi &; Mac Low! (120061 ) and iLi et al.l (120081 ) have studied the properties 



of tur bulence with ambipolar diffusion in a two-fluid approximation, while iPadoan et al. 



( I2000l ) have investigated the heating through ambipolar diffusion in turbulent molecular 
clouds using a single-fluid approximation. 

Shock waves in molecular clouds spread into steady-state continuous shocks, or C-type 
shocks, through ambipolar diffusion. If the shock speed is slower than the ion Alfven speed 
but faster than the neutral sound speed, the ions coupled to magnet ic fields drag the neu- 
trals into the postshock region, producing a continuous structure (jDraind Il980h . These 
steady-state C-type sh ocks can , howe ver, be unstable on a short enough timescale to be of 
astrophysical interest. IWardld (Il99ll ) showed that if the magnetic field lines are perturbed 
slightly and ions collect in the magnetic valleys, the ion-neutral friction may overcome the 
magnetic forces in the shock front and derive an exponentially growing instability. 

Numerical treatments of ambipolar diffusion have been commonly derived from ideal 
magnetohydrodynamic (MHD) models. Extensive numerical methods including ambipolar 



diffusion have been proposed in the stud y of the dynamics of partia 



the frame of single or two fluid models ( Toth 1 1 9 94 



19971 : ISmith fc Mac LoJl997l : lstone!ll997 



Li et al. 



Mac Low et al. 



2006: 



y ion i zed plasmas within 



1995 



Tilley fc Balsara 



Viae Low fc Smith 



20081 ). Mac Low et al 
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( 119951 ) have descr ibed an expl icit method for one-fluid ambipolar diffusion in the strong cou- 
pling limit, while iTothl (119941 ) has used a semi- i mplicit scheme for two- f luid ambipolar diffu- 
sion to investigate instability in C-type shocks. iTilley fc Balsaral (120081 ) have also presented 
a semi-implicit method for ambipolar diffusion using a two-fluid approximation. Implicit 
schemes for the multiflu id tr eatment of Hall diffus i on and amb ipolar diffusion have been 
suggested bv iFallel hoO$ ) and I O' Sullivan & Downesl (bood . l2007h . 



In this work we describe a fully explicit method for incorporating the single-fluid am- 
bipolar diffusion into a multidimensional MHD code based on the total variation diminishing 
scheme. The divergence-free condition of the magnetic field is ensured by a flux-interpolated 
constrained transport scheme, and a super time stepping method is used in order to consid- 
erably accelerate the otherwise painfully short diffusion-driven time steps. 

The organization of this paper is as follows. In §2 the MHD equations are presented 
along with the approximations we have made and in §3 our numerical methods are described 
in detail. Two test problems are presented in §4, while MHD turbulence simulations with 
significant ambipolar diffusion follow in §5. A summary is given in §6. 



2. MHD Equations with Ambipolar Diffusion 



We assume the strong coupling approximation, i.e., that the ion pressure and momentum 
are usually negligible in the weakly ionized plasma compared to those of the neutrals and 
so the magnetic force on the ions and the drag force exerted by the neutrals on the ions are 
almost equal. In this approximation the plasma can be represented as a single fluid. This 
single-fluid approximation turns out to be u seful in the form ation of molecular cloud cores 
through the process of ambipolar diffusion (IShu et al.l 119871 ). To simplify the modeling of 
ambipolar diffusion here we assume isothermality with a constant sound speed in the ions 
and the neutrals and we ignore gravity. 



The isothermal MHD equations including ambipolar diffusion can then be written as 

(1) 
(2) 



f + V- M = 0, 



dv 11 
— + v • Vv + -Vp (V x B) xB = 0, 

ot p p 



dB 

~dt 



Vx (v x B) = Vx 



1 



IPiP 



V x B) xB 



xB 



with the additional constraint for the absence of magnetic monopoles, 

V • B = 0. 



(3) 



(4) 
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Here the equation of state is p = a 2 p, where a is an isothermal sound speed, 7 is the collisional 
coupling constant between ions and neutrals, and pi is the ion density. The other variables 
p, v, and B denote neutral density, neutral velocity, and magnetic field, respectively. We 
renormalize the magnetic field by defining B = B/^/4tt throughout this paper so that the 
factor of 47r does not appear in equations (2) and (3). 

Although the ion density in molecular clouds depends on complicated physical balance 
between the cosmic-ray ionization of neutrals and the recombination of ions and electrons 
on charged grains, for the purpos e of simplicity we assume that the ion density scales as a 



power of the neutral density (e.g.. lElmegreenlll979l ) 



Pi = pio 



A) 



(5) 



By setting a = we further simplify by taking the ion density constant in this work, as these 
variations will not significantly affect the problems we treat. For conditions appropriate to 
molecular clouds, however, the choice of a ~ 0.5 usually would be more realistic, with 
a ~ most applicable at very high densities where grains are the main charge carriers 
(e.g., 



Ciolek &: Mouschoviaslll998l ). On the other hand, the ion velocity is obtained from the 



equation for the relative drift velocity between ions and neutrals, 



1 



Vi 



V 



ipip 



[V x B) xB. 



(6) 



This equation shows that the drag force and the magnetic force on the ions are balanced and 
that the ion-neutral drift velocity Vd = Vi — v is always perpendicular to the magnetic field. 

The basic effect o f ambipolar diff usion on the magnetic field can be expressed in a 
diffusion coefficient D (jShu et al.lll987l ) given by 



D 



where r = l/"fPi is the mean collisional time between ions and neutrals and ca 
the Alfven speed. Then we can estimate the ambipolar diffusion timescale as 



(7) 

B/y/p is 



D ' 



where L is the characteristic length scale of magnetic field. 
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Numerical Methods 



The numerical scheme 



works (IRvu et al.lll995l . Il998 



3.1. Source Term Integration 

or solving the id eal MHD equations is described in previous 



Kim et al. 19991 ). This method is based on the total variation 



diminishing (TVD) scheme (IHartenlll983l ) which is an explicit Eulerian upwind scheme with 
a second-order accuracy in space and time. In the strong coupling approximation we can 
separate the ion density and velocity from the neutral density and velocity so that we can 
basically use the MHD TVD code to compute the evolution of the neutral density and 
velocity using mass and momentum conservation equations. 

We now describe how to incorporate ambipolar diffusion terms into the induction equa- 
tion. The induction equation can be rewritten in component form as 



dB T d 



dy 



at 

dBy 
dt 

dt 



{B x v y 



ByV x ) 



d_ 

dz 



{B z v x - B x v z 



d d 

— (ByV Z - B Z Vy) - — (B X Vy ~ B yV , 



9 (n 



9 (R 

dy ^ 



B Z Vy) 



0S Z 


dSy 


dy 


dz ' 


dS x 


dS z 


dz 


dx ' 


dSy 


dS x 



dx 



dy 



(9) 
(10) 
(11) 



Here the source term components are given by 
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12) 
(13) 
(14) 

Note that the source terms on the right-hand side of equations (9) to (11) have the same 
divergence form as the flux components on the left-hand side. 

Standard second-order finite difference operators are applied to the explicit discretiza- 
tion of the source components. Here we define the source components at grid centers, S x ^j^ 
S y ,i,j t ki and S Zt ij t k, while the n-th components of the TVD flux vectors in each direction, 
/i^Vi/2,j,fc' fy\i,j+i/2,k' anc l /i!i,i,fc+i/2> are defined at face centers. The first four components 
of the TVD flux vectors, through are the upwind fluxes associated with the trans- 
port of mass and momentum, i.e., mass and momentum advection fluxes, and the last three 
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components of the TVD flux vectors, through f( 7 \ represent the components of the elec- 
tric field. While keeping second-order accuracy, the source components can then be included 
in the TVD flux components as follows 

+ g (S z ,i,j,k + S z ,i+lJ,k) j (15) 
— 2 (*~VM> fc "I" Sy,i+l,j,k) i (16) 



f (6) 


... f(6) 


f (7) 

J x,i+l/2,j,k 


_ f (7) 

J x,i+l/2,j,k 


f (7) 

Jy,i,j+l/2,k 


f (7) 

Jy,i,j+l/2,k 


AS) 

Jy,i,j+l/2,k 


Jy,i,j+l/2,k 


A5) 

J z,i,j,k+l/2 


_ f (5) 
J z,i,j, k+1/2 


f (6) 

J z,i,j,k+l/2 


Jz,i,j,k+l/2 



1 

2 



(17) 



(19) 
(20) 



2 \^Vthh k ~^ Sy,hh k + 1 ) ' 

2 \^x,i,j,k ^x,i,j,k+l) ■ 

Since the TVD scheme has second-order accuracy, the above second-order interpolation of 
the source components should be adequate. The contributions of the source components are 
added in a fully unsplit way after all the TVD flux components are updated through the 
TVD step. These total advective fluxes at face centers, /^ +1/2iJ>fc , fy% + i/ 2 ,k> and /i^,fc+i/2> 
are used to enforce V • B = as described in the following subsection as well as to update 
the magnetic field components to the next time step. 



3.2. Divergence-Free Condition 



Analytically the divergence- free condition V • B = is maintained if the condition holds 
for the initial magnetic field, but numerically the divergence of the magnetic field will not be 
exactly zero due to numerical discretization and dimensional splitting. Several sche mes to 
maintain V • B = cons t raint have been suggested and used in numerical MHD (see iToth 



20001) . Evans fc Hawleyl (119881 ) suggested the constrained transport (CT) scheme which 



used a specific finite difference discretization on a s taggered mesh to satisfy the divergence- 
free constraint. The flux-interpolated CT schemes (IRyu et al.lll998t iBalsara fc Spicerlll999l ) 
have introduced a new staggered magnetic field variable which is updated by simple finite 
differences using the interpolated fluxes. We have found the flux-interpolated CT scheme to 
be effective for incorporati ng ambipolar diffusion and have followed the approach suggested 
by lBalsara fc Spicerl (119991 ). 
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Using the total fluxes at face centers, equations (15) to (20), the advective fluxes at grid 
edges are reconstructed with second-order accuracy as follows 

O , , - - ( f ( ) I/I 6 ) _ f( 7 ) _ f( 7 ) ^ (ou 

"x,»j+l/2,fc+l/2 ^ \J z,i,j,k+l/2 ^ J z,i,j+l,k+l/2 Jy,i,j+l/2,k J y,i,j+l/2,k+l J ' 

o , , - I f f ( ? ) -i- f ( ? ) _ f ( 5 ) _ f ( 5 ) ^ roo^i 

"j/,j+l/2j,fc+l/2 — ^ ^J x ,i+l/2J,fe ^ J x,i+l/2 J,fe+1 Jz,i,j,k+l/2 J z,i+l,j,k+l/2j ' 

O , , - 1 f f ^ -i- f ( 5 ) _ f ( 6 ) _ f ( 6 ) ^ f9^ 

"z,i+l/2,j+l/2,fc ^ ^2/,ij+l/2,fe " r Jy,i+l,j+l/2,k ■> x,i+l/2j,k J x,i+l/2,j+l,k J ' v^ / 

Then the magnetic field components at face centers are updated as 

& "Jfl/2,j,fc = K,i+l/2,j,k ~ £- {^z,i+l/2,j+l/2,k ~ ttz,i+l/2j-l/2,k) 

+ ^ {^y,i+l/2,j,k+l/2 - ^y,i+l/2,j,k-l/2) , (24) 

^,ti+l/2,fc = b y,i,j+l/2,k ~ ^ fox,i,j+l/2,k+l/2 ~ ^x,i,j+l/2,k-l/2) 

+ ^ foz,i+l/2,j+l/2,k - ^2,j-l/2j+l/2,fe) , (25) 

& "ti,fc+l/2 = & "i,j,fc+l/2 - ^ {^y,i+l/2,j,k+l/2 ~ ^-1/2,^^+1/2) 

+ (^x,i,j+l/2,fc+l/2 - ^x,iJ-l/2,fe+l/2) • (26) 

It is straightforward to show that V • b n+1 = V • b n = if the numerical divergence of b is 
initially zero. 

In the TVD scheme for MHD, all fluid quantities are defined at grid centers. Thus the 
magnetic field components at grid centers are interpolated as 

Bx,i,j,k = 2 (px,i+i/2,j,k + b x ,i-i/2,j,k) , (27) 

B y,i,j,k = 2 ipy,i,i+i/2,k + &j/,i,j-i/2,fc) , (28) 

B z ,i,j,k = 2 (bz,i,j,k+i/2 + b z ,i,j,k-i/2) ■ (29) 

Note that the above arithmetic interpolation will be sufficient to maintain second-order 
accuracy. 



- 8- 



3.3. Super Time Stepping 



The time step for ambipolar diffusion is proportional to the square of the grid size in the 
single fluid appro ximation, so the explic it treatment of ambipolar diffusion terms leads to very 
small time steps (IMac Low et al.lll995l ). In the two fluid approximation that treats ions and 
neutrals separately, including the ion momentum equation severely limits the time steps via 
a very restrictive stability criterio n. To resolve t his problem, previous papers have proposed 
different solutions. For instanc e. iLi et al.l (120061 ) prop osed a "heavy-ion" approximation to 
speed up the time steps, while iNakamura fc Lil (120081 ) set a density threshold below which 
the ambipolar diffusion rate is s et to zero to avoid thi s problem. In this work we adopt the 
"super time stepping" approach (jAlexiades et al.lll996f ) to incre ase the effective time i n terval 



and allow much faster computations for ambipolar diffusion. |Q' Sullivan Downed (12006 
20071 ) also used this strategy in their multifluid MHD models. 



The super time stepping techn i que c onsiderably accelerates the explicit schemes for 
parabolic problems (jAlexiades et al.lll996l ). The key advantage of this approach is that it 
demands stability over large compound time steps, rather than over each of the constituent 
substeps. In this method, the state vector is evolved over a super time step, 



At 



STS 



At 



n 



AD 



2v^ 



l + v^) 2n -(l- >fr) 



in 



1 



v^) 2 " + (1 - v^) 

where A£ad is the nominal ambipolar diffusion timestep, n is the number of substeps, and 
v is a fuzzy factor (0 < v < 1). This super time step is defined as Atgxs 
consists of the substeps, Atj, which are given by 

-l 



>2n 



(30) 



J21=i Atj, and 



Ai 



At 



AD 



{v-l) 



COS 



2j-ln 



n 



+ v + l 



(31) 



It has been proven that AtsTS _ > ^ 2 A£ad as v — > so that the su per time step approac h 



is asymptotically n times faster than the standard explicit scheme (jAlexiades et al 



19961 ). 



However, the v parameter must be properly chosen for each problem in order to achieve 
optimality and stability of performance. For diffusion-dominated problems, A£ad drops 
below the typical Courant time step, At C our (i-e., At ad < Ai S TS < At Cour ). Thus, if At ST g 
is taken to be the Courant time step, the super time stepping method requires roughly 
A/At Cour /At A D substeps. 

In addition to allowing larger effective time steps, the super time stepping approach 
offers relatively simple implementation since it is a first order method. We have successfully 
applied this approach to the linear and nonlinear ambipolar diffusion problems presented 
in the following sections. Our numerical results confirm the efficiency and accuracy of the 
super time stepping approach, as previously implemented in other approximations. 
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4. Test Problems 

4.1. Oblique C-type Shocks 

We compute the structure of oblique C-type shocks in order to test the numerical meth- 
ods described in the previous section. The steady state structure of a C-type shock is 
characterized conveniently by the shock length scale defined as 

L s = rc A . (32) 



Following iMac Low et al.l (119951 ). the steady state equations are solved by setting d t — d y — 
d z = in equations (1) to (3), which then can be reduced to a single ordinary differential 
equation for p. The steady state solution can be obtained through the numerical integration 
of the ordinary differential equation. In a separate code we numerically integrate the ordinary 
differential equation using the fourth-order Runge-Kutta method. This solution is specified 
by the three parameters, the sonic Mach number M, the Alfven Mach number Ma, and the 
angle 9 between the shock normal and the magnetic field. 

To generate C-type shocks we set up a two-dimensional shock heating problem. Initially 
a gas with a neutral density p propagates with a velocity v against a reflecting wall placed 
at x = in the uniform magnetic field B that lies at an angle 9 to the x-axis. As the 
gas hits the reflecting wall, both the fluid and the magnetic field are compressed, the gas 
is heated and a reverse shock is produced. The ion-neutral friction drags the neutral gas 
into the postshock region, and finally the steady-state C-type shock is built up, yielding the 
appropriate continuous transition. The parameters we chose for this problem are a = 0.1, 
7 = 1, pi = 10~ 5 , p = 1, and B = B 0x x + B 0y y with B 0x = B 0y = l/\/2. This problem has 
been set up with two inflow velocities, v = —4A5x and — 9.47^, which correspond to the 
shock velocities v s = 5 and 10 for our chosen parameters. This gives M = v s /a = 50 and 
100, Ma = v s /cx = 5 and 10, and 6 = tan _1 (i?oj / /-Boa:) = 7r/4. The computations have been 
done in a two-dimensional box of x — y — [0, L] with L = 20L S using 128 2 cells. Outflow 
boundary conditions are used except for a reflecting boundary imposed at x — 0. 

The parameters for the run shown in Figure [T](a) are M = 50, Ma = 5, and 9 = tt/A and 
those for the run shown in Figure H^b) are M = 100, Ma = 10, 9 = tt/4. In Figured! the 
structure of neutral density, neutral (red) and ion (blue) velocity components, and magnetic 
field components from numerical calculations are marked with open circles and compared 
to analytic solutions plotted with solid lines. Structures are measured along the x-direction 
before the shock reaches the outer boundary. The spike in the neutral density seen in a 
few cells near the reflecting wall is the overheating phenomenon. This is purely a numerical 
artifact that most finite difference schemes applied to the shock heating problem intrinsically 
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possess since they cannot compute the jump condition across strong shocks within a single 
cell. In all the flow variables, the structure of the C-type shock clearly forms. Figure [1] shows 
the excellent agreement between the numerical solutions and the analytic solutions for the 
steady-state C-type shocks, demonstrating the accuracy of our numerical methods. 

The accuracy of numerical solutions depends on the number of cells spanned by the box 
size L. So we have run the case of the test in Figure Ufa) with different numerical resolutions 
to check the convergence properties. Except for the resolutions the initial conditions are 
identical to those used in the test in Figure Ufa). We have computed the mean errors 
for neutral density defined by E(p) = Ylij\Pi,j ~ Pi,j\/Tlij\Pi,j\i where the superscript 
n represents numerical solution and the superscript a represents analytic solution. The 
resolutions of 16 2 , 32 2 , 64 2 , 128 2 , and 256 2 cells give the mean errors of 1.33, 1.29, 1.28, 
1.27, and 1.27, respectively. As expected, the mean errors for neutral density converge as 
the numerical resolution increases. In this convergence test we also see that there are clear 
trends toward convergence in mean errors for velocity and magnetic field, and that a 128 2 
grid is sufficient to treat this problem. 



4.2. Decay of Alfven Waves 



The propagation of Alfven waves in a we akly ionized plasma provi des an effective tool for 
testing the dynamics of ambipolar diffusion. iKulsrud fc Pearcd ( 119691 ) first showed that am- 
bipolar diffusion can prevent the prop agation o f Alfve n waves in a partially ionized medium. 
In the strong coupling approximation iBalsaral (119961 ) gives an explicit quadratic dispersion 
relation for Alfven waves, 



on 2 + 

Ipi 



c a h 



0. 



(33) 



where to = ur + iuji is the complex angular frequency of the wave and k is a real wavenumber. 
It is clear in the above equation that the Alfven waves always propagate when k < 2'jpi/cA 
(i.e., ujr 7^ 0). In order to test the propagation of Alfven waves in the strong coupling 
limit, we have followed the evolution of a standing wave in numerical calculations and com- 
pared the oscillation frequency and decay rate of the wave t o the analytic results. Damped 
oscillations of standing waves have been long studied (e.g., Morse fc Ingardlll986l ) and the 
time-dependence of the first normal mode is described by 



h (t) = h | sin (ujRt) \ e' 



(34) 



where ho is the initial amplitude of the wave. 



In our test of the decay of Alfven waves, we have used a standing wave formed along 
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the diagonal on x — y plane with initial velocity 

v = v amp c A sin (k x x + k y y) z. (35) 

The background density and magnetic field have been set to be uniform with p = 1 and 
B = B x with B = 1. This gives the characteristic Alfven speed ca = B/y/2p = 0.7071. 
Here the initial peak amplitude has been set to v amp = 0.1 and the wavenumbers have been 
set to k x = k y = 2tt/L so the total wavenumber is k = ^Jk% + k 2 . We choose collisional 
coupling constants, 7 = 100, 500, and 1000 with a = 1 and pi = 0.1 in this test. The 
calculations have been done in a computational box of x — y — z — [0, L] with L = 1 using 
128 3 cells. Boundary conditions are periodic in the x- and y-directions and outflow in the 
^-direction. 

Figure [2] shows the time evolution of the (spatially) root mean square magnetic field in 
the ^-direction, < bB\ > 1//2 , for three different collisional coupling constants 7 = 1000 (top), 
500 (middle), and 100 (bottom) in the test of the decay of Alfven waves. Our numerical 
results are marked with open circles while the theoretical predictions from equation (34) are 
plotted as solid lines. The oscillation frequencies and decay rates from these numerical calcu- 
lations fit very well to those from the theoretical predictions, confirming that our numerical 
methods are accurate. Based on the time evolution of the standing wave in Figure El the 
propagation of Alfven waves is significantly suppressed with decreasing collisional coupling 
constant 7. 

By fitting the time evolution of < bB\ > x l 2 to theoretical curves from equation (34) we 
find the numerical data for the complex angular frequency u, whose real and imaginary parts 
correspond to the oscillation frequency and decay rate of the standing wave, respectively. We 
repeated the calculation of Figure[2]for 7 = 100, for eight wavenumbers ranging from 2tt / 16 to 
6V2n, collecting the data for the oscillation frequencies ujr and the decay rates uji. In Figure 
E]we show those "experimental" results together with the analytic solution of equation (33). 
Oscillation frequencies (red) and decay rates (blue) found from the numerical experiments 
are represented with filled circles and the analytic solutions are drawn as solid lines. In 
Figure [3] the very good agreement between numerical data and theoretical predictions for 
different wavenumbers shows the accuracy and flexibility of our numerical methods. 



5. MHD Turbulence Simulations 



In this section we present, as a first practical problem using this code, simulations of 
the compressible MHD turbulence in the pres ence of ambip o lar di ffusion. A simulation of 
turbulent ambipolar diffusion was studied by iPadoan et al.l (120001 ) in the strong coupling 
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approximation and by lOishi fc Mac Lowl (120061 ) and iLi et al.l (120081 ) in the two-fluid ap- 
proximation. We have performed this simulation to confirm the validation of our numerical 
methods for including ambipolar diffusion as well as to investigate the role of ambipolar 
diffusion in the dissipation of compressible MHD turbulence. 

To characterize the MHD simulations, we define the strength of the magnetic field in 
terms of the parameter 

a 2 



(36) 



Note that our definition of (3 differs by a factor of 2 from the usual plasma (3, the ratio 
of gas pressure to magnetic pressure. For MHD turbulent flows with the root mean square 
velocity v rms , the importance of magnetic fields on the dynamics of gas is characterized by the 
Alfven Mach number Ma = t>rms/ c A = \/~/3M, where the sonic Mach number M is given by 
M = v vms /a. The effect of ambipolar diffusion on A lfvenic turbule nce of scale L is measured 
through the ambipolar diffusion Reynolds number (jBalsaralll996l ) defined by 



R 



AD 



(37) 



A sufficiently large ambipolar diffusion Reynolds number implies that the importance of 
ambipolar diffusion to the turbulent flow on the scale L becomes vanishingly small. The 
ambipolar diffusion length scale /ad is then defined as a characteristic length scale at which 
the ambipolar diffusion Reynolds number becomes unity, i.e., /ad = i~c\/v rms . 

We con sider two types of MHD turbulence models driven according to the method 



described in IStone et al.l (119981 ). One is turbulence decaying from saturated initial velocity 
perturbations, and the other is forced turbulence in which velocity perturbations are added at 
constant time intervals. In both decaying and forced turbulence, the velocity perturbations 



5v are generated from a Gaussian random field with a power spectrum 

P (k) = 5v 2 (k) oc k 6 exp (Sk/k p ) , 



(3* 



where the power spectrum peaks at k p = 4(2ir/L). The velocity perturbations are subject to 
the constraints that V-<5u = and no net momentum is added by the velocity perturbations, 
< p5v >= 0. The perturbations are normalized so that the initial kinetic energy 5Ek = 50 
for decaying turbulence and a constant kinetic energy input rate E K = 500 is injected at 
regular time intervals 5t = 0.001 for forced turbulence. 

The simulation parameters for both decaying models and forced models are summarized 
in Table [TJ According to the values of 7 varying from 00 to 100, we denote decaying models as 
Dl to D3 and forced models as Fl to F3. In Tabled], the flow time is defined as tf = L/v rms , 
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and tend is the simulation end time. The model simulations have been set up with uniform 
neutral density p — 1 and uniform magnetic field B = B x with B — 1. The isothermal 
sound speed and ion density are assumed to be constant with a = 1 and pi = 0.01. The 
simulations have been done in a periodic box of x = y = z = [0, L] with L = 1 using 128 3 
cells. 

In Figures H](a) and (b) we present images of the logarithms of the neutral density for 
models Dl (top) and D3 (bottom) at t = It/ and for models Fl (top) and F3 (bottom) at 
t = ltf, respectively. The images are slices through the x — y plane at z = 0.5. The overall 
density features for both decaying turbulence and forced turbulence models are roughly 
the same, even though the turbulence driving pattern is different. Small-scale knots and 
filaments are produced in the absence of ambipolar diffusion (models Dl and Fl). By 
introducing strong ambipolar diffusion (i.e., reducing the collisional coupling constant to 
10), the density structures diffuse out, resulting in larger, smoother density structures, as 
shown in models D3 and F3. The results for models D2 and F2 are intermediate and are 
not shown. We note that the ranges of observed densities are smaller in the cases where 
ambipolar diffusion is strong, regardless of whether the turbulence is decaying or forced. 

Figure [5] shows the time evolution of the total energy, defined as the sum of the kinetic 
energy and the energy of the perturbed magnetic field, E tot = Ek + Eb, where 

1 r o. 



Ek = ^ / pv dV, (39) 



v 
and 

Eb = \J v {B 2 - Bl) dV, (40) 

with dV = dxdydz. The evolution of the total energy for the decaying models, Dl (black), 
D2 (blue), and D3 (red) are plotted in Figure 0(a). After initial plateau phases all these 
models lose their initial energies rapidly, decaying nearly as power-laws with time (with 
indices between 1.1 to 1.3). The rate of the turbulent energy decay significantly increases 
with decreases in the collisional coupling constant, and thus we see that the turbulent decay 
rate can be strongly affected by differences in ambipolar diffusion. The somewha t more rapid 



declin e of the total energy found for model Dl than for a similar simulation by IStone et al. 



(119981 ) can be understood in terms of the slightly lower resolution and higher initial energy 
in our simulation. In Figure [5](b), the evolution of the total energy for the forced models, 
Fl (black), F2 (blue), and F3 (red) are shown. In all the forced models the total energy 
rises steeply and then saturates at final states since the dissipation rate balances the input 
power. The amplitude of the final saturated energy level decreases with decreases in the 
collisional coupling constant, again showing that the dissipation rate increases as the strength 
of ambipolar diffusion increases. 
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6. Summary 

In this paper we describe specific numerical methods for incorporating ambipolar dif- 
fusion into a multidimensional MHD code based on the total variation diminishing scheme. 
We assume the strong coupling approximation, that magnetic force and the neutral-ion drag 
force in weakly ionized plasmas are almost equal and so the plasma can be treated as a single 
fluid. Since our numerical methods described in this paper are fully explicit and maintain 
a second-order accuracy, it is straightforward to extend them to parallelized versions and 
other geometries. The divergence-free constraint on the magnetic field has been exactly en- 
forced through the flux-interpolated constraint transport scheme at all times. By using the 
super time stepping method to accelerate the timestep for ambipolar diffusion, we remove 
the severe restriction on the stable timestep that would arise at high numerical resolution 
and/or in strong ambipolar diffusion. 

Ambipolar diffusion has been tested through the direct comparison with analytic solu- 
tions of diffusion problems. We have computed test problems that include oblique C-type 
shocks and the decay of Alfven waves. For both of these test problems, comparisons of 
numerical results to analytic solutions are possible and they demonstrate the good accuracy 
and robustness of our methods. We have also performed simulations of the compressible 
MHD turbulence in the presence of ambipolar diffusion and they confirm the ability of our 
code to follow complex MHD flows. We have shown that the dissipation rate of MHD turbu- 
lence is strongly affected by the strength of ambipolar diffusion in both decaying turbulence 
and forced turbulence. 

This multidimensional MHD code incorporating an explicit scheme for solving the am- 
bipolar diffusion term allows us to study astrophysical systems such as molecular cloud cores 
and protostellar discs in which ambipolar diffusion is thought to be important. Currently 
this code is being used to study the evolution of compressible MHD turbulence with am- 
bipolar diffusion, and results of three-dimensional, high resolution MHD simulations will be 
reported elsewhere. 

This work utilized a high performance cluster at the Korea Astronomy and Space Science 
Institute (KASI). EC was supported by the KASI Postdoctoral Fellowship. JK was supported 
by the Korea Science and Engineering Foundation through the Astrophysical Research Center 
for the Structure and Evolution of Cosmos and by the Korea Foundation for International 
Cooperation of Science and Technology through grant K20702020016-07E0200-01610. 



- 15 - 



REFERENCES 

Alexiades, V., Amiez, G., k Gremaud, P.-A. 1996, Comm. Num. Meth. Eng., 12, 31 
Balsara, D. S. 1996, ApJ, 465, 775 

Balsara, D. S., & Spicer, D. S. 1999, J. Comput. Phys., 149, 270 

Basu, S., & Ciolek, G. E. 2004, ApJ, 607, L39 

Ciolek, G. E., & Mouschovias, T. Ch. 1998, ApJ, 504, 280 

Desch, S. J., & Mouschovias, T. Ch. 2001, ApJ, 550, 314 

Draine, B. T. 1980, ApJ, 241, 1021 

Elmegreen, B. G. 1979, ApJ, 232, 729 

Evans, C. R., & Hawley, J. F. 1988, ApJ, 332, 659 

Falle, S. A. E. G. 2003, MNRAS, 344, 1210 

Fatuzzo, M., & Adams, F. C. 2002, ApJ, 570, 210 

Harten, A. 1983, J. Comput. Phys., 49, 357 

Heitsch, F., Zweibel, E. C, Slyz, A. D., & Devriendt, J. E. G. 2004, ApJ, 603, 165 

Kim, J., Ryu, D., Jones, T. W., & Hong, S. S. 1999, ApJ, 514, 506 

Kulsrud, R., & Pearce, W. P. 1969, ApJ, 156, 445 

Li, P. S., McKee, C. F., & Klein, R. I. 2006, ApJ, 653, 1280 

Li, P. S., McKee, C. F., Klein, R. L, & Fisher, R. T. 2008, ApJ, 684, 380 

Mac Low, M.-M., Norman, M. L., Konigl, A., & Wardle, M. 1995, ApJ, 442, 726 

Mac Low, M.-M., & Smith, M. D. 1997, ApJ, 491, 596 

Mestel, L., & Spitzer, L., Jr. 1956, MNRAS, 116, 503 

Morse, P. M., & Ingard, K. U. 1986, Theoretical Acoustics (Princeton: Princeton Univ. 
Press) 

Mouschovias, T. Ch. 1976, ApJ, 207, 141 



-16- 

Mouschovias, T. Ch., & Ciolek, G. E. 1999, in The Origin of Stars and Planetary Systems, 
ed. C. J. Lada & N. D. Kylafis (Dordrecht: Kluwer), 305 

Myers, P. C. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews (Tucson: 
Univ. Arizona Press), 81 

Nakamura, F., & Li, Z.-Y. 2008, ApJ, 687, 354 

Oishi, J. S., & Mac Low, M.-M. 2006, ApJ, 638, 281 

O'Sullivan, S., & Downes, T. P. 2006, MNRAS, 366, 1329 

O'Sullivan, S., & Downes, T. P. 2007, MNRAS, 376, 1648 

Padoan, P., Zweibel, E., & Nordlund, A. 2000, ApJ, 540, 332 

Ryu, D., Jones, T. W., & Frank, A. 1995, ApJ, 452, 785 

Ryu, D., Miniati, F., Jones, T. W., & Frank, A. 1998, ApJ, 509, 244 

Shu, F. H., Adams, F. C, & Lizano, S. 1987, ARA&A, 25, 23 

Smith, M. D., & Mac Low, M.-M. 1997, A&A, 326, 801 

Stone, J. M. 1997, ApJ, 487, 271 

Stone, J. M., Ostriker, E. C, & Gammie, C. F. 1998, ApJ, 508, L99 

Tassis, K., & Mouschovias, T. Ch. 2007, ApJ, 660, 370 

Tilley, D. A., & Balsara, D. S. 2008, MNRAS, 389, 1058 

Toth, G. 1994, ApJ, 425, 171 

Toth, G. 2000, J. Comput. Phys., 161, 605 

Wardle, M. 1991, MNRAS, 251, 119 

Zweibel, E. G. 2002, ApJ, 567, 962 



This preprint was prepared with the AAS IATgX macros v5.2. 



-17- 




Fig. 1. — (a) Structure of an oblique C-type shock with M = 50, Ma = 5 and 9 = n /4. 
Profiles of neutral density, neutral (red) and ion (blue) velocity components, and magnetic 
field components from the numerical calculation are marked with open circles. The solid 
lines represent the analytic solution for the steady-state C-type shock. The calculation has 
been done in a square box with size L = 20L S using 128 2 cells, (b) Same as in (a) except for 
M = 100, M A = 10, and 9 = vr/4. 
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Fig. 2. — Time evolution of the spatially averaged root mean square magnetic field in the in- 
direction for 7 = 1000 (top), 500 (middle), and 100 (bottom) in the test of the decay of Alfven 
waves. Our numerical results (open circles) are compared to the theoretical predictions (solid 
lines). The calculations have been done in a cube box with size L = 1 using 128 3 cells. Time 
is expressed in units of the sound wave crossing time, L/a. 
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Fig. 3. — Oscillation frequencies ur and decay rates ujj collected from the oscillation of 
standing waves with different wavenumbers for 7 = 100 in the test of the decay of Alfven 
waves. Oscillation frequencies (red) and decay rates (blue) are marked with filled circles and 
the analytic solutions are drawn as solid lines. The calculations were done in a cubic box of 
size L = 1 using 128 3 cells. 
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Fig. 4. — (a) Images of the neutral density for models Dl (top) and D3 (bottom) at t — ltf 
in the MHD turbulence simulations. The simulations have been done in a cubic, periodic 
box with size L = 1 using 128 3 cells, and the images shown are slices through the x — y 
plane at z = 0.5. The color bars are drawn in logarithmic (base 10) scales, (b) Same as in 
(a) except for models Fl (top) and F3 (bottom) at t — ltf. 
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Fig. 5. — (a) Time evolution of the total energy for the decaying models Dl (black), D2 
(blue), and D3 (red) in the MHD turbulence simulations. The simulations have been done 
in a cubic, periodic box with size L = 1 using 128 3 cells. Time is expressed in units of the 
sound wave crossing time, L/a. (b) Same as in (a) except for the forced models Fl (black), 
F2 (blue), and F3 (red). 
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Table 1. Parameters for MHD Turbulence Simulations 



Model 


P 


7 


M A 


T 


-Rad 


^AD 


*/ 


*end 


D1(F1) 


1 


oo 


10(1) 





oo 





0.1(1) 


10(1)*/ 


D2(F2) 


1 


1000 


10(1) 


0.1 


100(10) 


0.01(0.1) 


0.1(1) 


10(1)*/ 


D3(F3) 


1 


100 


10(1) 


1 


10(1) 


0.1(1) 


0.1(1) 


10(1)*/ 



Note. - Here models Dl to D3 represent decaying turbulence and 
models Fl to F3 denote forced turbulence. All the decaying and forced 
models have been done in pairs, and the different values for both models 
are distinguished using parentheses. 



